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Abstract 

This paper addresses a multi-scale finite element method for second order linear elliptic equations 
with arbitrarily rough coefficient. We propose a local oversampling method to construct basis func¬ 
tions that have optimal local approximation property. Our methodology is based on the compactness 
of the solution operator restricted on local regions of the spatial domain, and does not depend on any 
scale-separation or periodicity assumption of the coefficient. We focus on a special type of basis func¬ 
tions that are harmonic on each element and have optimal approximation property. We first reduce 
our problem to approximating the trace of the solution space on each edge of the underlying mesh, 
and then achieve this goal through the singular value decomposition of an oversampling operator. 
Rigorous error estimates can be obtained through thresholding in constructing the basis functions. 
Numerical results for several problems with multiple spatial scales and high contrast inclusions are 
presented to demonstrate the compactness of the local solution space and the capacity of our method 
in identifying and exploiting this compact structure to achieve computational savings. 

Keywords. Multi-scale finite element method. Oversampling. Optimal local basis. High-contrast. 


1 Introduction 

Many problems of practical importance in science and engineering have multi-scale feature: composite 
materials modeling and flows in porous media are typical examples of such kind. In many cases, the 
qualities of interest are only related to the large-scale properties of the solutions. However, the fine-scale 
features of the model can have significant impact on the large-scale properties of the solutions, thus 
one needs to use very fine mesh to resolve the small-scale variations of the problem even though only 
large-scale properties of the solutions are required. The computational cost can be prohibitive. For these 
so-called multi-scale problems, efficient upscaling methods that allow us to incorporate the small-scale 
features of the problem into the large-scale properties of the solutions are desired. 

In this work, we use the following example of second order linear elliptic equation with homogeneous 
Dirichlet boundary condition to illustrate our upscaling methodology, 

f —div(a(*)Vu(a:)) = f(x), x G SI; 

\u(aO|an = 0, 

where fl is a convex polygon domain in R d with d = 2,3. We assume that the equation is uniformly 
elliptic, i.e., there exist A m j n > 0 and A m i n > 0 such that 

fl(x) G [A m in, A max ]. (1-2) 

We do not assume any regularity of the coefficient a(x) G L°°(ff), which may have multiple spatial scales, 
thus the above equation (1.1) can be used to model diffusion process in strongly heterogeneous media. 
We also assume that in (1.1) the forcing function f(x) G L 2 (Q). Then the existence of solution to (1.1), 
u(x ) G follows immediately from the Lax-Milgram theorem, and we have 

cll/IL-qn) < Hz)lltfi(n) < Cll/lln-ipi)- (1.3) 

However, due to the roughness of the coefficient a(x), the solution to (1.1) u(x ) loses regularity and 
ceases to be in H 2 (Cl). Classical finite element method uses piecewise polynomials to approximate the 
solution space, and its convergence depends on the following approximation property 

IML - Ju(*)||ifi ( n) < Ch\\u(x)\\ H ^ n) , (1.4) 

* hou@cms.caltech.edu, plliu@cms.caltech.edu. Computing-(-Mathematical Sciences, Caltech. 


1 



and the regularity result 


IWz)lli/ 2 (fi) < C'||/(a:)|| L 2 (fl) , (1.5) 

where Ju is the piecewise polynomial interpolation of u(x), and h is the underlying mesh size. Classical 
finite element method may fail for these problems (1.1), since ||u(a:)||jj 2 cannot be bounded by ||y(;c)||z^ 2 (rj) 
in (1.5). It is actually showed in [6] that the polynomial finite elements can perform arbitrarily badly in 
this setting, thus (1.1) can serve as a typical example of multi-scale problems. 

One of the strategies to numerically solve the multi-scale problem (1.1) when classical finite element 
methods fail is using problem-dependent basis that incorporates properties of the coefficient a(x), to 
approximate the solution space. To be specific, one constructs basis functions 

01 {x),(j> 2 (x),...<l>n(x) £ (1.6) 

that may depend on the elliptic operator —div(a(x)V(-)), and find numerical solution 

u h {x) £ Vh(x) = span{0i(x), 0 2 (x),... 0 n (x)} C Hj(fi), (1.7) 

using the Galerkin projection. Namely, we find Uh(x ) within the trial space Vh, such that 

a(u h (x),v(x)) = (f(x),v(x)), for all v £ V h , (1.8) 

where 

a(u(x), v(x)) = / Vu(x) t a(x)Vv(x)dx, (/( x),v(x)) = / f(x)v(x)dx. (1.9) 

Jn J n 

The numerical solution defined above satisfies the following optimal property under the energy norm 

\\u(x) - u h {x)\\ a = inf ||u(a;) — v(x)\\ a , (1.10) 

v(x)ev h 

where the energy norm is equivalent to the Hq(Q) norm, and defined as 

!l u ( a: )lla = a (u(x), u(x)) = / Vu(x) t a(x)'Vu(x)dx. (1.11) 

J n 

In this work we will employ the above strategy to numerically solve (1.1). Note that to obtain the 
numerical solution Uh{x) from the Galerkin projection (1.8), one needs to solve a linear system of size 
n x n. Thus to make the computational cost small, we want the number of the basis functions used in 
(1.6) to be small. Besides, we want the basis functions in (1.6) to have compact support such that the 
stiffness matrix formed in (1.8) is sparse thus easy to compute and invert. 

We propose an effective method to construct basis functions (1.6) with optimal local approximation 
property in this paper. Our method is based on the compactness of the solution operator to (1.1) 
restricted on local regions of the domain. To be specific, we introduce the following operator 

Ti : f(x) -» Ui(x) = u(x)\ Di , (1.12) 

where D t is local subset of fl of size 0(H ), and H is chosen according to the desired order of accuracy. 
We construct local basis functions that can approximate the range of Tj with controlled accuracy, and 
combine them together to get the approximate solution space to (1.1). The compactness of the operator 
Ti, after removing the singular part, will be demonstrated numerically in section 2. 

On each local region of the domain, Di, we divide the local solution Ui(x) into two orthogonal parts 
with respect to energy norm (1.9): an o(x)-harmonic part, and a local bubble part. We show that the 
bubble part of the solution is small and its compact structure can be easily identified by inverting the 
elliptic equation (1.1) locally on each region Di. We consider approximating the a(a;)-harmonic part 
of the solution space using a special type of basis functions that are a(a;)-harmonic on each Di (but 
not across the boundary of Di), and call basis functions of such type multi-scale basis. Due to the 
smallness of the bubble part of the solution, we demonstrate that multi-scale basis functions are optimal 
in approximating the solution space for fixed local boundary conditions on dDi. 

The u(x)-harmonic part of the solution only depends on the restriction of the solution on the boundary 
of the local regions Di, and we seek to identify the compact structure of the trace of the solution space 
on dDi. Using a primary set of interpolation basis functions, i/>i(x), we can reduce our problem to 
approximating the solution space on each edge of the coarse mesh, e. We introduce an oversampling 
operator that maps the solution on an oversampling domain W to the solution restricted on an edge of 
the mesh e. Then we employ the compactness of the oversampling operator to construct the optimal 
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boundary basis functions on each edge e. The optimal choice of the interpolation basis functions ipi(x) 
can also be identified by solving least square problems. With these local boundary basis functions, we 
construct basis functions (1.6) that approximate the a(x)-harmonic part of the solution space by solving 
some local boundary value problems. The resulting basis functions (1.6) are a(a:)-orthogonal with respect 
to the bubble part of the solution space, and because of this property we can add the bubble part back 
to our numerical solutions by simply solving some local cell problems. 

The resulting method consists of two stages: in the offline stage we identify the local compact structure 
of the solution space, and build multi-scale basis functions and the corresponding stiffness matrix in (1.8); 
in the online stage, for any given forcing function f(x) £ L 2 (fl), we solve the equation (1.1) efficiently 
using the multi-scale basis functions constructed offline with a very low computation cost. Our method 
can achieve significant computational savings in the multi-query setting where equation (1.1) need to be 
solved for multiple times with different forcing functions f(x) £ A 2 (12). 

Numerical examples for several problems with rough coefficients and high-contrast channels are pre¬ 
sented. Our method achieves high accuracy for these problems, and these numerical results suggest that 
our method is very robust and works equally well for problems without scale separation, or have high 
contrast channels. We demonstrate our methodology through the second-order scalar elliptic equation, 
but it can be applied to other linear elliptic problems like elasticity equations without much difficulty. 

There is a fast growing literature on multi-scale methods for elliptic equations with rough coefficient, 
and we briefly review some related work below. The classical homogeneziation theory [7] considers 
a(x) with periodic structure a(x) = a(x, -), and derives an effective equation governing the asymptotic 
behaviors of the solutions as the small scale e goes to 0. In [21, 18, 22, 14], the authors proposed the 
Multi-scale finite element method (MsFEM), in which the multi-scale basis functions are constructed by 
solving local elliptic problems with homogeneous right hand side, and the convergence analysis of MsFEM 
in the periodic setting was given in [18, 15]. An oversampling technique [18] was proposed to reduce the 
resonance error introduced due to the artificial boundary conditions in solving the local problems. This 
present work can be viewed as a continuation of the MsFEM, and we propose a robust method to choose 
and enrich the local boundary conditions using the oversampling operator. In [27], the authors showed 
that the solutions gain an order of regularity with respect to the Harmonic coordinate [3], and constructed 
multi-scale basis functions using this property. The Harmonic coordinate was recently employed in [11] for 
global upscaling. In [8, 28], the flux norm was introduced and employed to show the compactness of the 
solution space and construct (localized) basis functions. In [29] the polyharmonic spline was employed, 
which was later put in the Bayesian inference setting [26]. In [24], the generalized finite element method 
was proposed, which provides a general framework to combine local approximation spaces together using 
a partition of unity formulation. In [5, 4], the local special basis functions are constructed by solving 
some local spectral problems, and then combined together using the partition of unity framework. The 
generalized multi-scale finite element method [12] is a systematic method to construct multi-scale basis 
functions for a family multi-scale problems with parameters. The Heterogeneous multi-scale method 
[1, 30] numerically decomposes the structure of the medium into a micro-scale and a macro-scale, and 
compute solutions of cell problems on the micro-scale to get the local homogenized matrices. Finally we 
remake that multi-scale methods using problem-dependent basis functions were also employed to solve 
problems in which the coefficient has high contrast inclusions [9]. 

The remaining part of this paper is organized as follows. In section 2, we demonstrate the compactness 
of the solution operator restricted on local regions of the spatial domain. In section 3, we divide the 
solutions on each local region of the domain to different parts corresponding to the trace of the solution 
on the edges of the coarse mesh, and identify their compact structures separately. And with this compact 
structure, we construct multi-scale basis functions that have optimal local approximation property. In 
section 4, numerical results are presented to demonstrate the capacity of our method in identifying and 
exploiting the compactness of the solution space to achieve computational savings. Concluding remarks 
are made in section 5. 


2 Compactness of the Solution Space Restricted on Local Re¬ 
gions of the Domain 

The existence of finite number of basis functions (1.6) that can approximate the solutions space to (1.1) 
up to any accuracy is implied by the compactness of the solution operator, T, which maps from the 
forcing function f(x) £ A 2 (12) to the corresponding solution u(x) £ A7q(12). 

T: f(x) e L 2 {Q) ^u{x) eH^{Q). (2.1) 
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The compactness of T is demonstrated in [23, 8], and it was employed for elliptic equations with random 
input data recently in [19, 20] for stochastic model reduction. 

To be specific, the solution operator T can be decomposed as 

T = L (2-2) 

where L maps f(x) £ H~ l (Q.) to the solution u(x) £ Hq( fi), and 1 (n) is the embedding 

operator from L 2 (Q) to ,ff _1 (fi). From (1.3), we can see that L is continuous and indeed a homo¬ 
morphism, and the compactness of is well known based on the Sobolev space theory [16]. 

Thus the compactness of T follows from the decomposition (2.2). To quantify the approximability of T 
by a finite-rank operator, we consider its Kolmogorov-n width, which is defined in below. 

Definition 2.1 (Kolmogorov n-width). For a compact linear operator T that maps between two Hilbert 
spaces, we define its Kolmogorov n-width as 


d n (T)=mf||r-T n ||, (2.3) 

where T n runs over all rank-n linear operators. 

Due to the fact that L is a homomorphism, one can easily see that the Komogorov-n width of T 
is only different from that of t(n) by a factor that depends on A m i n , A max (1.2) and D, 

cd n (I L ?( f2)_>,H--i(fi)) < d n (T) < Cd n (I L 2(n)_>/f-i(o)). (2.4) 

The Kolmogorov-n width of the embedding operator is well-known [23, 8], and we have 

d n (lL 2 (ci)->H^(a)) = n 1 ^ d (C + o( 1)), n —* oo. (2.5) 

From (2.4) (2.5) and (2.3), we obtain that there exist of a set of basis functions, (1.6), with the following 
approximation property to the solution space of (1.1), 

n 

sup inf II 'y'c i <t>i{x) - w(x)|| H i (f2) < Cn~ 1/d . (2.6) 

ll/ODMnj-a c< ~i 

The approximation property (2.6) is optimal, and does not depend on the regularity of the coefficient 
a(x). For practical applications in multi-scale problems, we want the basis functions 4>i(x) to have local 
support such that the corresponding stiffness matrix in (1.8) is sparse and easy to invert. However, the 
basis functions in (2.6) whose existence is implied by (2.4), (2.5) and (2.6) may be nonlocal. 

Since our objective is finding basis functions (1.6) with local support, we consider a local region of 
the domain, D with diameter O(H), and a slightly larger local domain which contains D, W, which we 
call the oversampling region. We consider the restriction of the solutions to (1.1) on W, 



Uw{x) = u(x)\w- 

(2.7) 

The local solution uw{x ) 

can be divided into two parts, 



u w {x) = Uwfa) + Uw{x), 

(2.8) 

where 

(— div(a(a:)Vw^(a:)) = 0, x £ W, 

|«^(a;) = u w [x), x £ dW, 

(2.9) 

and 

f—div(«(s:)Vt4(i))=/(j:), x £ W, 

|it^(x) = 0, x £ dW. 

(2.10) 


We call the first part u\y(x) the local a(x)-harmonic part, and the second part u^y(x) the local bubble 
part. The two parts are orthogonal with respect the local inner product, aw(', ')i 


/ \/u\ v {x) t a{x)^Jul v {x)dx = 0. (2.11) 

Jw 

The local bubble part u^y(x) is small in the sense that 

\\ u w ( x )\\ 2 h £( w ) — < ^'-^“ll/( a: )ll! 2 (W’)’ (2-12) 
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which can be obtained from (1.3) and a scaling argument. (2.12) implies that if we only want to obtain 
O(H) accuracy in our numerical solution, we can simply neglect the local bubble part. 

Then we consider a local solution operator Tp that maps f(x) to the local a(a:)-harmonic part, uly(x) 
restricted on D, 

Tp : f(x) £ L 2 (Q) i4 (*)|a e H\D), (2.13) 

and we want to construct local basis functions on D that can approximate the range of Tp. To demon¬ 
strate the compactness of To, we choose a set of orthonormal basis in the domain and range of To to 
discretize Tp as a matrix, and compute the decay of its singular values. We consider the following choice 
of coefficient in (1.1), which has multiple fine spatial scales and is illustrated in Figure la, 


a(x) 


1 

6 


1.1 + sin(27rr/e 1 ) 

1.1 + sin(27r?//e 1 ) ~*~ 


1.1 + sin(27rj//e 2 ) ( 1.1 + cos(27nr/e 3 ) 

1.1 + cos(2-kx / e?) 1.1 + sin(27n//e 3 ) 

1.1 + sin(2irj//e 4 ) 1.1 + cos(27nE/e 5 ) 
1.1 + cos(27nr/e 4 ) 1.1 + sin(27rj//e 5 ) 


+ sin(4ari/ 2 ) + 1 


(2.14) 


where e, = §, e 2 = e 3 = vf, £4 = e 5 = ^ 

We choose fi = [0,1] x [0,1], the oversampling region W = [1477, 1777] x [1477, 1777], and the local 
region D = [1577,1677] x [1577,1677], where 77 = 1/32. The decay of the singular values of the local 
solution operator (2.13) is plotted in Figure lb. Then we compute the singular values for the local 
solution operator (2.13) to the Possion equation in the same setting, the decay of which is plotted in 
Figure lc. From the figures, we can see that the singular values of the local solution operator decay very 
fast, and this fast decay does not deteriorate due to the roughness of the coefficient. 



(a) The rough coefficient. 


(b) Low rank, (rough coefficient) 


(c) Low rank. (Possion equation) 


Figure 1 


The fast decay of singular values of To implies that we can use a very small number of local basis 
functions, to be specific, the first several left singular vectors of Tp, to get very good local approximation 
property. However, we cannot afford to construct Tp explicitly since it is a solution operator and its 
construction involves solving the equation (1.1) a large amount of times (globally). It is known that for a 
low-rank operator, the main action oiTp can be captured in its image on some random vectors. This, to 
some degree, explains the success of some global upscaling methods [13, 11, 27] that use sampled global 
solutions to (1.1) to approximate the solution space to (1.1) locally. 

We will not pursue this perspective in this work. Instead, we decompose Tp using a global operator 
and a local oversampling operator, and construct local multi-scale basis functions employing the com¬ 
pactness of the oversampling operator. The resulting method does not involve any global solving of the 
equation (1.1), and will be detailed in the next section. 


3 Identify the compact structure of the solution space using 
Oversampling 

In this section, we identify the compact structure of the local solution space through oversampling, and 
use it to construct basis functions (1.6). In our numerical examples, the domain ff is chosen to be 
[0,1] x [0,1], and we discretize S 2. using a coarse square mesh of size 77, which should be chosen according 
to the desired order of accuracy. With this discetization, we have 

fi = u£r D u (3.1) 
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where Di have disjoint interiors. Underlying this coarse mesh, we use a triangle fine mesh of size O(h), 
which is a refinement of the coarse mesh. The fine mesh size h should be chosen such that it can resolve 
the small scale variation of the multi-scale coefficient in (1.1). In our method we solve the equation (1.1) 
on the coarse mesh, and the basis functions that we use are constructed and saved using linear basis 
functions on the fine mesh. The two level discretization is illustrated in Figure 2. 

Two Level Mesh 



Figure 2 


3.1 The multi-scale basis 

We first introduce a special class of problem-dependent basis functions (1.6), which we call the multi¬ 
scale basis. 

Definition 3.1 (Multi-Scale basis). For a discretization of 11 (3.1), we consider basis functions 


4>i{ x )i <t>2(z), ■ ■ ■ <t> n {x) e -ffo(fi). (3.2) 

If they are a(z)-harmonic on each coarse element of the coarse discretizaion, Dj , 

— div(a(z)V<)>i(a;)) = 0, x G Dj, (3-3) 

then we call them multi-scale basis functions. 

Clearly, the multi-scale basis functions are determined by their restrictions on the boundary of coarse 
elements dDi since they are a(z)-harmonic in each Di. We denote 

T = U t^Di- (3.4) 

We have the following proposition, which implies that if the desired accuracy is O(H), multi-scale basis 
functions are optimal for fixed local boundary conditions on T (3.4). 

Proposition 3.1. Consider a set of basis functions ipifx) G HQ(fl),i = 1,2,. ..m, and a set of multi¬ 
scale basis functions (j>i{ x) G Hq (fl), i = l,2...,n <ra u coarse mesh of size H, as showed in Figure 2. 
Denote the corresponding Galerkin numerical solution (1.8) to (1.1) using tfi(x),i = 1 ... .m as Uh{x), 
and the Galerkin solution using <j>i(x), i = 1,.. .n as uff s (x). If 

spaniel(z)| r , • . .0»(a:)|r, ■. .^„(x)| r } = span{V>i(*)|r, • ■■tfi{x) |r ■. . ip m {x)\r}- (3.5) 

Then vj€- hcivG 

ll«(z) - Uh S (x)\\l < Hz) - u h (x)\\l + C||/||| 2(f7) IF 2 . (3.6) 

Namely, if only O(H) accuracy in the energy norm is desired, the multi-scale basis can perform as well 
as other set of basis functions, given that the local boundary conditions of basis functions are the same. 





























To prove the above proposition, we first divide the solution u(x) to (1.1) into two parts. On each 
coarse mesh element Di , we consider 

Ui(x) = u(x)\ Dt , (3.7) 

and divide it to an a(x)-harmonic part and a local bubble part, as we did in (2.8), 

Ui(x) = u}(x) + uf(x), x e Di, (3.8) 

where u\(x) is the local a(a)-harmonic part, and u 2 (x) is the local bubble part. Combine these local 
decompositions from all coarse elements Di together we get, 

N N 

u(x) = u}(x) + U 2 (x), v}(x) = Uj(x), U 2 (x) = u i( x ). (3-9) 

i =1 i= 1 

One can see that the two parts v}(x), u 2 (x) are orthogonal with respect to the a (•, •) inner product (1.9), 

a(u 1 (x), u 2 (x)) = 0. (3.10) 

Besides, the combination of the local bubble parts is small according to (2.12). To be specific, 

\\u 2 (x)\\ a <CH\\f\\ L2{n) . (3.11) 


Next we prove the proposition 3.1. 

Proof. Denote the numerical solution using ipi(x),i = 1 ,...m, Uh(x) as Uh(x) = YlrLi^i^iix), tlien, 
there exist Ci,i = 1,.. .n such that u™ s = X^=i <k<l>i(x), an d 

ur(x)\r = u h (x)\ T . (3.12) 

Then we consider 

\\Uh 3 ( x ) ~ u ( x )\\l = ll« 2 (z) + m 1 ^) - «“ s (*)llo- (3T3) 

Since u™ s (x) £ .ffg(fi) is o(a:)-harmonic on each coarse element Di, we have 

a(u 2 (x),u™ s (x)) = Y'' f '77u 2 (x) t a(x)'77u™ s (x)&x = — f u 2 (x)div(a(x)'Vu]( l3 (x)dx = 0. (3.14) 

i.,i >», 

Thus u 2 (x) is a-orthogonal to u 1 (a;) — u™ s (x), and according to (3.11) we have 

H*) - < S (x)\\l = \\A x )fa + h\x) - ur(x)\\l < llu^ar) - n^x)\\ 2 a + C\\f\\ 2 LHn) H 2 . (3.15) 

Then we consider u e (x) = u(x) — Uh(x ), and divide it to two parts as we did for u(x ) in (3.9). We get 


u e (x) = u\(x) + u 2 (x), a(ug(x), u 2 (x)) = 0. (3.16) 

Consequently, we have 

IK*) - u h (x)\\l = \\ul(x)\\l + \\ul(x)\\l > \\ul(x)\\l. (3.17) 

According to (3.12), we have 

u\(x) = u l (x) - u^ix), (3.18) 

since they are equal on T and a(x)-harmonic on each Di. 

Finally based on (3.15), (3.17), and the optimal property (1.10), we have 

ll«(*) - u h S (x)\\l < IN*) - < s {x)\\l < \W(x) ~ u h (x)\\ 2 a + C\\f\\ 2 L 2 (n) H 2 , (3.19) 

and complete the proof. □ 


As we have showed in (3.11), the bubble part of the solution u 2 (x) is small and of O(H) in the 
energy norm, thus can be neglected if the desired accuracy in the numerical solution is O(H). In our 
method in this work, we use multi-scale basis functions in (1.6) to approximate the solution space, which 
are locally a(x)-harmonic functions, and are a(x)-orthogonal to the bubble part of solution. Due to 
this a(a:)-orthogonality and the Galerkin projection formulation in (1.8), multi-scale basis functions only 
approximate the o(*)-harmonic part of the solution and will not bring in additional errors in the bubble 
part. Thus we can recover the bubble part of solution u 2 (x) independently by solve some local bubble 
problems (2.10). And adding u 2 (x) back to uff s (x), we can get numerical solution that is free of error 
in the bubble part. This is one of the advantages of using multi-scale basis functions in (1.6). 

To construct local multi-scale basis functions, we divide the o(x)-harmonic part of the solution v}(x) 
into different parts corresponding to different edges of the coarse mesh, and approximate them separately. 
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3.2 Decomposition of the a(x)-harmonic part of the solution 


To identify the compact structure of the a(a;)-harmonic part of the solution, we first introduce a set of pri¬ 
mary interpolation multi-scale basis ij)i(x), i = 1,. .. n based on the coarse mesh node points X\, * 2 , • • • x n , 

ipi(xj) = 5ij\ — div(a(x)'Vipi(x)) = 0, x £ Dj. (3.20) 

We also require that ij)i(x) is supported on the four coarse elements around aFor example, we can 
simply choose the multi-scale basis ipi(x) to be linear on the boundaries of coarse elements. We will 
discuss about the optimal choice of these primary interpolation basis functions in subsection 3.4. 

For f(x) G L 2 (f2), and the spatial dimension d = 2, 3, we have that u(x) is Holder continuous on Q 
[17], so we can consider the interpolation of u 1 (*), namely the a(x)-harmonic part of the solution, using 
the primary basis functions and get the residual, 

ul(x) = u 1 ^)-^2u(xi)ipi(x). (3.21) 

i 

For a coarse mesh element D t , we denote its four nodes points as x j 15 Xi 2 , x» 3 and ®j 4 , then we get 
the restriction of the residual (3.21) on Di , 

u\{x)\ Di = u](x) - uixD^ix) - u(x i2 )ij)i 2 (x) - u(x i3 )ip i3 (x) - u(x u )ip u (x). (3.22) 


n 





_ e.i _ Xi 2 


e 2| 

A b 


Xi 3 . . 

_£4 _ I^i 4 


(a) Decomposition of the interpolation residual. 


Oversampling Region W 


X 


£1--£3-fc-ft)- X 




D u 


Target edge e 

D d 

-£»- '*4 


-£ie- 


^5 


PT 


(b) Oversampling region. 


Since the residual Mg(x)|o. vanishes on the node points Xi j: j = 1,2,3,4, we can divide the trace of 
u] on dDi to four parts, corresponding to the four edges of H,, ei, e 2 , e 3 , e±, 

«e(*) \d Di =u\{x)\ ei + u\ {x)\ e2 + ul(x)\ e3 + u\{x)\^. (3.23) 

This decomposition is illustrated in Figure 3a. Each part in the above decomposition (3.23) belongs to 
H l f 2 {dDi) since they vanish on the node points, and we can extend them to Di to get four a(x)-harmonic 
components of ul(x). We denote them as v ei (x), v e2 (x), v e3 (x), v ei (x), and have 

ul(x)\ dDi = V ei (x) + Ve 2 (x) + Ve 3 (x) + V e 4 (x). (3.24) 

Combining these local decompositions together, we have 

ul(x) = J^v e (x), (3.25) 

e 

where v e (x) is the o(x)-harmonic extension of the interpolation error on the edge e to its two neighbor 
elements. In the above decomposition (3.25), we are actually dividing the error uf (x) in the a(x)-harmonic 
part of solution into different parts corresponding to errors on different edges e. This is possible since 
uf(x) vanishes on the node points, thanks to the interpolation operation using ipi(x) (3.21). 

We seek to construct boundary basis functions on each edge e that approximate v e (x), and combine 
them together to get the whole trial space. We introduce the following operator for the edge e with 
endpoints x 4l and Xi 2 , which maps f(x) G L 2 (fl) to the interpolation residual of the solution on e, 

T e : f(x) G i 2 (fi) -» v e (x) = u(x) - u(x il )(l>i 1 (x) - u(xi 2 )<j>i 2 (x) G H 1/2 (e). (3.26) 

The left singular vectors of T e form the optimal local boundary basis functions. However, T e is a 
global operator and its construction involves solving the equation (1.1) globally. In the next section, we 
decompose T e as a global solution operator and a local oversampling operator, and construct boundary 
basis functions that approximate the range of T e through the oversampling operator. 












3.3 The oversampling operator 


To identify the compact structure of the solution space restricted on the edge e, we put it in an oversam¬ 
pling region that we denote by W. In our numerical examples, we use the square mesh of size H for the 
coarse discretization, and the oversampling region W is chosen as the union of the six elements around 
the edge e. It is illustrated in Figure 3b. We remark that our method is also applicable to other types 
of discretizations like triangular mesh. We denote the solution on W as Uw(x). 

We remark that the idea of identifying the local structure of the solution space by putting it in 
a larger region, namely oversampling, was first proposed in [18] to reduce the resonance error due to 
artificial local boundary conditions, and this strategy was later employed in [2, 5, 10]. 

We denote T\y as the operator that maps f(x) £ L 2 (Q) to the oversampling solution uw{x) = u(x)\w, 
and T\v~ye as the operator that maps uw(x) to the solution restricted on the edge e: 

T w : f(x)^u w (x) = u(x)\ w , T w ^ e : u w (x) -» u e (x) = u w (x)\ e . (3.27) 

We also introduce the interpolation residual operator using (3.20), P e , 

P e : u{x)\ e ^u{x)\ e -u(x il )i> il (x)-u(x i2 )ip i2 (x). (3.28) 

With the above definitions, the operator T e (3.26) can be decomposed as 

T e = P e T w ^ e T w . (3.29) 

We call the operator P e Tw^e hr the above decomposition (3.29) the oversampling operator, which 
maps the solution on W , uw{x) to the interpolation residual, which we denote by v e (x), 

Pos = PeTw^e ■ u w (x) -s- v e {x) = u w (x) - u w (x il )ip h (x) - u w (x i2 )i> i2 (x), (3.30) 

where x 2l and *j 2 are the two endpoints of e, and ipi(x) is the primary interpolation basis (3.20). 

The solution to (1.1) on W, uw(x ) can be divided into two parts, the a(a:)-harmonic part u\y(x) and 
the bubble part u^ v (x). We employ the compactness of the oversampling operator (3.30) to construct 
basis functions in H 1 ^ 2 (e) that vanish at x 2l and Xi 2 , and approximate the range of (3.26). To be specific, 
we use the first several left singular vectors of Pos as the basis functions associated with e. We first 
introduce appropriate inner products for the domain and range space of Pos- 

On the edge e, the image of T e , v e {x) £ H 1 ^ 2 (e) and vanishes on the two endpoints. We consider its 
a(*)-harmonic extension to the upper and lower coarse elements respectively, as showed in Figure 3b, 
and denote them as v™(x) and v^(x). Then we define 

\\v e (x )\\ 2 H1/ 2 (e) = \ J D Vv e (x) l a(x)Vv"(x)dx + V v*(xf a(x)V v^(x) da:. (3.31) 

In the domain of the operator P e Tw^ei namely, uw{x), we define its inner product as 

\\ u w(x)\\l= [ Vu] v (x) t a(x)\/u] v (x) + (w^(x)) 2 dx + [ [di v(a(x)\7u w (x))] 2 . (3.32) 

Jw Jw 

With the above inner products, we compute the singular value decomposition of the oversampling 
operator Pos- To discretize the domain of Pos, we consider its two parts, the a(x)-harmonic part, and 
the bubble part. The a(a;)-harmonic part only depends on the trace of u\\r{x) on dW, we discretize 
ff 1 / 2 (dW) using all the fine mesh piecewise linear functions. If dW intersects with dQ, then we choose 
H 1 ^ 2 (dW) to be fine mesh basis functions that vanish on <9S2. The bubble part of the solution u^ v (x) 
only depends on fw(x) = f(x)\w, and we discretize f(x) using piecewise constant functions on the 
coarse mesh. The following lemma justifies this discretization of f(x) £ L 2 (Q). 

Lemma 3.1. Denoe the space of piecewise constant functions on the coarse mesh as V c . The coarse 
mesh is showed in Figure 2 with mesh size H. Then for f(x) £ L 2 (Q), there exist f c (x) £ V c , such that 

||/( 2 :) - / c (x)||if-i ( n) < CH\\f{x)\\ L 2 (n) . (3.33) 

According to (1.3), the above Lemma implies that if we want to obtain O(H) accuracy in the energy 
norm of the numerical solution, we can only consider piecewise constant forcing functions on the coarse 
mesh. If higher accuracy is desired in the numerical solution for f(x) £ L 2 (fl), we can simply refine the 


9 



discretization of L 2 (fl) in the oversampling operator (3.30). We remark that if f(x) has higher regularity, 
we can get better approximation result in (3.33), for example 


inf 

/cGV c 


!!/(*) - fc(x)\\ H -i {n} < CH 2 \\f(x)\\ H i (n) . 


(3.34) 


With the above discretization of u\v(x), we truncate the singular values of Pos to O(H) and select 
the corresponding left singular vectors as the boundary basis functions e. We denote them as 

vl(x),...v*‘(x) e f? 1/2 (e), (3.35) 

which vanish on the two endpoints of e. According to the oversampling operator in (3.29), and our 
truncation criteria, we have the following approximation property 

efc 

inf ||u(a:)|e - ufajfaix) - u(x i 2 )<j> i 2 (x) - ^c>*(x) \\ H u^( e ) < CH( ||uN*)IL + \\f\\ L *(w))- (3-36) 

‘ i= 1 


Then we extend these boundary basis functions to the two neighbourhood coarse elements D u and 
D d as a(a;)-harmonic functions, by solving 


—div(a(x)'V4>e{x)) =0,x e D u ,D d , 


(3.37) 


~ u e V 


Finally, we combine the multi-scale basis functions associated with each edge of the coarse mesh, e, 

<j>l(x), i = 1,2... fc e . (3.38) 

with the primary interpolation basis functions ipiix) (3-20), and get the trial space, 

V h = span{<(>j(a:), i=l,...n}. (3.39) 


We have the following error estimate using the trial space (3.39). 

Proposition 3.2. Using the trial space consisting of the primary interpolation basis functions (3.20) and 
the multi-scale basis functions constructed from each edge of the coarse mesh (3.37), we obtain numerical 
solution to (1.1) using the Galerkin projection (1.8). Then we have the following convergence property, 

IN*) - wf S (*)lU < CH\\f(x)\\ L 2 . (3.40) 


Remark 3.1. Using a simple Aubin-Nitsche duality argument and the convergence result in the energy 
norm (3.40), we can get the following convergence result in the L 2 norm. 

\\u(x) - ut IS (x)\\ L 2 (n) < CH 2 \\f(x)\\ L 2 . (3.41) 


The proof of the convergence result (3.40) follows directly from the decomposition of the solution 
operator (3.29) and the truncation in the singular value decomposition of the oversampling operator. 


Proof. We choose cj, as the ones in (3.36), and denote 

n k e 

< s {x) = y, u(Xi)lpi(x) + EE ^(x) e V h . (3.42) 

i =1 e j =1 

Then we consider \\u™ s (x) — u(a;)|| a , since the basis functions in (3.39) are multi-scale basis, namely, they 
are o(a)-harmonic on each Di , we have 

ll«r(*) - u(*)lla = ll«r(*) - Mx) - n 2 (x)\\l < IICW - «i(*)lla + CH 2 \\f\\l Hn) , (3.43) 

where u\(x) and u 2 (x ) are the a(a:)-harmonic part and bubble part of the solution u(x). 

Then we divide — ui(x)||^ to different parts on Di. For each part, according to the approxi¬ 

mation property (3.36), and the definition (3.31), we have 



v(«r(*)-M*))N*)v(«r(*) 


u i(*))d* < CY H \\ u w(x)\\l- 

w 


(3.44) 
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The sum over W corresponds to the oversampling regions for edges of D 2 . There are four of them for 
each Di. Summing up (3.44) for all the coarse elements of Q, we have 

\\< s (x) ~ MaOHo ^ CH 2 \\u(x)\\ a . (3.45) 

Putting (3.45) in (3.43), we have 

ll«(*) - CWIU < CH\\f(x)\\v> {a) . (3.46) 

Then using the optimal approximation property (1.10), we finish the proof. □ 

To make the number of multi-scale basis functions in (1.6) small, we want the singular values of 
Pos = PeTw^e decay fast. Pos can be divided into two parts: the first part acts on the a(a;)-harmonic 
part of uw(x), and we denote it as Pqs'j the second part acts on the bubble that only depends on fw{x), 
and we denote it as Pqs- For the first part, similar analysis has been done in [5] in a slightly different 
setting, and the method there also applies to our problem. We have the following result. 

Proposition 3.3. Denote the singular values of Pq S that acts on a(x)-harmonic functions on W as <r k , 
then we have the following upper bound on the decay of a k : for any e > 0, there exist C such that 

a k < Cexpl—fc 1 /( d+1 )~ e }. j (3.47) 


where d is the dimension of the domain SI. 

The second part Pq S is small according to (2.12), and the decay rate of its singular values can be 
obtained from (2.5) using a simple scaling argument. 

We will see in our numerical results section that the singular values of Pos decay very fast, and a 
very small number of boundary basis functions can achieve high local approximation accuracy. 

As we have shown previously, the basis functions we obtain are multi-scale basis functions, thus are 
a(s)-orthogonal to the bubble part of the solution space. This gives us the flexibility to add the bubble 
parts back to the numerical solutions at local regions where higher accuracy is desired by simply solving 
some local bubble problems (2.10). In our truncation of the singular values of the local compact operator 
PeTw^te i we choose the threshold to be O(H), since O(H) accuracy is required in (3.40). If we need 
higher accuracy than H, for example, O(e) where 

h « e « H. (3.48) 

Then we can truncate the singular values of P e Tw^ e by e, by doing which, the resulting multi-scale basis 
functions are able to approximate the o(x)-harmonic part of the solution space up to 0(e) accuracy. Then 
by adding back the bubble part of the solution u 2 (x) to the numerical solution uff s (x), we can get O(e) 
accuracy in our final numerical solutions. Namely, our upscaling strategy allows us to get arbitrarily 
high accuracy that is permitted by the fine mesh discretization. 


3.4 Optimal primary interpolation basis functions 

In constructing the multi-scale basis functions in the previous section, we need to choose a set of primary 
interpolation basis functions (3.20) first, which allows us to reduce the problem to approximating the 
solutions space restricted on each edge e, (3.29). The choice of these interpolation basis will affect the 
oversampling operator Pos = PeTw-te- In this subsection, we identify the optimal interpolation basis 
functions by solving local under-determined least square problems. 

The oversampling operator for edge e, Pos , depends on the interpolation basis functions tpi 1 (x) and 
ipi 2 (x) associated with the two endpoints of e, and we seek optimal primary interpolation basis functions 
and ipi 2 (x) that make the singular values of Pos have the fastest decay. 

We consider the following optimization problem 

min HV’u(a:)||a> min \\i/)i 2 (x)\\ a , subject to (3.49a) 

-div(a(x)VV»i 3 .(*)) £ L 2 {W), x 6 W, ipij{x ik ) = S jk - (3.49b) 

where the norm j| • || a is defined in (3.32). 

Under certain conditions, the optimization problem (3.49) has unique solutions, and the unique 
solutions are optimal interpolation basis functions. We have the following theorem. 
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Theorem 3.1. Let(j>* l ( x),(j>* 2 {x) be the solution to (3.49), and(j>i 1 (x) and (j>i 2 (x) be two other multi-scale 
interpolation basis functions. Denote the corresponding oversampling operators using these interpolation 
basis as (3.30) as Pq S and Post an d let <r*,i = 1,2,... and a.i,i = 1,2,... be their singular values. 
Assume that both a* and o, are arranged in descending order. Then we have 

<ri < (3.50) 

which implies the interpolation basis in (3.49) make the singular value of the oversampling operator (3.30) 
have the fastest decay. 

The proof of Theorem 3.1 is given in the appendix. 

Note that in the minimization problem (3.49), the two interpolation basis functions can be constructed 
independently by solving under-determined least square problems. By going over the oversample regions 
for each edge of the coarse mesh, we can construct the optimal interpolation basis functions (3.20) on all 
the boundaries of local regions. Then we can extend them locally to o(a:)-harmonic functions by solving 
some local boundary value problems as in (3.37) to get the interpolation basis (3.20). 

We will see in our numerical results section that, in some cases, the optimal interpolation basis 
functions (3.49) themselves are enough to approximate the solution space of (1.1). Namely, there is no 
need to construct the multi-scale basis functions associated with each edge (3.37). 

We also remark that the optimization problem (3.49) can be interpreted as a Bayesian Inference 
problem as in [29, 26, 25]. Interested readers may consult the references for more information. 


3.5 Implementation of the whole method 

The proposed multi-scale finite element method consists of two stages, the offline stage and the online 
stage. In the offline stage, we identify the compact structure of the solution space. In the online stage, 
for a given forcing function, we compute the numerical solution using the offline basis functions. 

The offline stage involves the following procedures, 

1. Build the a(x)-harmonic extension operator on each oversampling region W. 

On each oversampling region IV, we build the local a(*)-harmonic extension operator, which maps 
the boundary condition which belongs to H 1 ^ 2 (dW) to a(a:)-harmonic functions on W. This step 
requires solving a series of boundary value problems on W. 

2. Construct the H 1 / 2 (dW ) norm. 

The a(ai)-harmonic extension operator and the inner product (3.32) together induces a norm on 
H 1 ^ 2 (W). We need this weighted norm in the SVD of the oversampling operator. 

3. Construct the optimal multi-scale interpolation basis functions, ipi ± (x) and ipi 2 (x). 

Discretize the domain of the oversampling operator Pos using the strategy described in section 3.3, 
and build the restriction operator that maps %(i) to WiyOUi) and uw(xi 2 ). Then we solve the 
optimization problem (3.49), which is essentially an under-determined least square problem. 

Solving the optimization problem (3.49) on each oversampling region only gives us the boundary 
conditions of the primary interpolation basis functions, we need to solve local boundary value 
problems to get their local o(x)-harmonic extensions as the interpolation basis functions. 

4. Build the oversampling operator. 

Using the optimal interpolation basis functions we constructed in the previous step, we build the 
oversampling operator Pqs , (3.3). 

5. Construct the H 1 / 2 {e) norm. 

Construct the H x l 2 (e) norm (3.31), which requires solving a series of local boundary value problems 
on the two neighbour coarse elements of e . 

6 . Compute the SVD of the oversampling operators. 

Using the inner product (3.32) and (3.31) in the oversampling operator to compute its singular 
value decomposition. We truncate the singular values to e, and save the corresponding left singular 
vectors, which are basis functions of i7 1 / 2 (e), vl(x), v 2 (x),... Vg e (x). 


12 



7. Construct multi-scale basis functions associated with each edge e. 

We extend the basis functions of H 1 ^ 2 (e), namely u*(x), to the two neighbourhood coarse elements 
of e by solving local boundary value problems (3.37). 

8 . Compute the stiffness matrix. 

Combining the multi-scale interpolation basis functions with the multi-scale basis functions asso¬ 
ciated with each edge of the coarse mesh, we get the trial space 

V h =spa,n{<j) 1 (x),<j>2(x),...<l> n {x)}. (3.51) 

Then we save these multi-scale basis functions, and compute the stiffness matrix, 

= a(4>i{x),4>j{x)). (3.52) 

The online stage involves the following procedures, 

1. Compute the load vector. 

For a given forcing function f(x) 6 L 2 (Q), we compute the corresponding load vector 

b(i) = f 4>i(x)f(x)dx, i = l,...n. (3.53) 

Jn 

2. Compute the online numerical solution. 

Using the load vector b (3.53) and stiffness matrix (3.52), we solve the linear system 

Me = b, (3.54) 

and get the online numerical solution on the fine mesh 

n 

Uh S {x) = £ciMx). (3.55) 

i=l 

Recall that the multi-scale basis functions associated with the edges vanish on the coarse grid node 
points. So if we only want the large scale solution, we can simply select the coefficients in c that 
correspond to the primary multi-scale interpolation basis (3.20). 

3. Recover the bubble part of the solution. 

Solve the local boundary value problem (2.10) on each coarse mesh element D {, and get u 2 (x). 
Combine these local bubble parts together and add them back to Galerkin solution (3.55), we get 

N 

Uh(x) = uff s (x) + ^ «?(*)• (3-56) 

i=l 

If 0(H ) accuracy is required in the numerical solution, this step is unnecessary since the bubble 
part does not impact the large scale properties of the solution. 

Note that in the offline stage, we need to solve a series of boundary value problems for each edge 
of the coarse mesh to construct the oversampling operator (3.30), and then compute its singular value 
decomposition, which is computationally expensive. However, the constructions of multi-scale basis 
functions on each edge are independent from each other, thus the offline stage can be implemented on a 
parallel machine to accelerate the computation. In the online stage, the main computational cost comes 
from solving the linear system (3.54). Our numerical results in the next section suggest that a small 
number of basis functions are enough to obtain the coarse mesh accuracy, O(H), thus the linear system 
(3.54) is small and sparse. This implies that the online computation in our method is efficient, and our 
method can bring in significant computational savings in the multi-query setting, where the equation 
(1.1) needs to be solved for multiple times using different forcing functions. 


4 Numerical Results 

In this section, we present numerical examples that have multiple-scale features and high-contrast chan¬ 
nels to demonstrate the capacity of our method in identifying and exploiting the compact structure of 
the local solution space to achieve computational savings (in the online stage). We discretize the domain 
of the problems fi = [0,1] x [0,1] using a two level mesh as showed in Figure 2. The coarse mesh is of 
size H = 1/32, and the fine mesh is of size h = 1/1024. 
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4.1 An example with multiple spatial scales 


The first example we consider is one that has multiple spatial scales. The coefficient is given by (2.14), 
and it is visualized in Figure la. For each edge of the coarse mesh, we compute the singular value 
decomposition of the oversampling operator, and truncate the singular values at e = H, which guarantees 
O(H) accuracy in the online numerical solution. After the offline stage, multi-scale basis functions (1.6) 
are constructed, and the average number of basis functions associated with each edge is 


ke 


E e fe e 

#(e) 


1 . 00 . 


(4.1) 


k e is very small. Actually only 1 or 2 multi-scale basis functions are constructed for each edge of the 
coarse mesh, aside from the primary interpolation basis functions. And this reflects the effciency of our 
method in the online stage since the stiffness matrix is small and sparse. 

To measure the error in our online numerical solution, we need to choose a reference solution. Since 
the multi-scale basis functions are constructed and saved on the fine mesh of size h, we will use the 
piecewise linear finite element solution on the fine mesh as the reference. 

In the online stage, we choose the forcing function f(x,y) to be 


f{x,y) = l, (x,y)eQ. 


(4.2) 


Recall that the basis functions that we use are a(a;)-harmonic in each Di , and we can add back the 
bubble part of the solution by simply solving some local cell problems. Our online numerical solution, 
and the numerical solution after correction using the bubble part are plotted in Figure 4. We can see 
that the numerical errors in the online numerical solution is very small. 



(a) The fine mesh solution. 


(b) Error in the numerical solution, (c) Error corrected using the bubble. 


Figure 4: Online numerical solutions 


We measure the error of the numerical solution in the energy norm and L 2 norm. We denote the 
numerical solution as ujj S (x), and the corrected solution using the bubble part as uh{x). We compute 


i MS 

\\u(x) -U% S (x)\\ a 

E„ = 

| «(z) - U H (x)\\a 

(4.3a) 

a 

\Hx)\\a 

!K*)IL 

iMS 

\\u{x) - u% s {x)\ |i,2 (n) 

El 2 = 

l|u(®) - U H (x)\\ L 2 ( n) 

(4.3b) 

' L 2 

IKz)IU 2 (f1) 

||u(aOI|z,2 ( fi) 


The results are listed in Table 1. We can see that by adding the bubble part of solution back to the 
numerical solution, the error in L 2 norm and energy norm are both reduced by about one half. This 
implies the numerical error in the o(a;)-harmonic part of the solution is about the same as that in the 
bubble part. The latter is of order O(H) in the energy norm and of order 0(H 2 ) in the L 2 norm, and 
this result confirms our error estimates (3.40) and (3.41). 

Then we consider only using the primary interpolation basis functions ifii ( x ), i = 1,... N in the trial 
space (3.51), namely, we do not use the multi-scale basis functions associated with each edge of the 
coarse mesh. The error in the corresponding numerical solution is also listed in Table 1. We can that 
the relative error in L 2 is also small, which means the numerical solution can capture the large-scale 
property of the solution. However, the errors in the energy norm and L 2 norm are both significantly 
larger than that in (4.3), which implies the necessity of enriching the trial space using the multi-scale 
basis functions associated with the edges of the coarse mesh if higher accuracy is required. 
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Energy Norm 

L 2 Norm 

Numerical solution 

4.16 x 10“ 2 

1.73 x 10~ a 

Corrected numerical solution 

2.67 x 10“ 2 

8.75 x 10~ 4 

Interpolation basis only 

7T5'7~x~TcP^ 

5.95 x 10- a 


Table 1: Numerical errors 


To further demonstrate the convergence rate in (3.40) and (3.41), we consider a sequence of coarse 
mesh size with H = 2~ k ,k = 3,4, 5,6,7. For each H , we compute the error in the online numerical 
solution, the decay of the numerical error with H is plotted in Figure 5. From the plot, we can clearly 
see that the decay rates of the online numerical error agree with the estimates (3.40) and (3.41). 




(a) Energy norm. 


(b) L 2 norm. 


Figure 5: Convergence of the online numerical solution with N = 1/H. 


4.2 An example without scale-separation 

In this subsection, we consider an example of the coefficient a(x) without scale-separation. 

a(x,y) = |a| +0.5, (4.4) 

where a is normally distributed on the node point of a mesh of size And a(x,y) on the fine mesh 
is obtained using piecewise linear interpolation. This coefficient a(x, y) is rough and has no clear scale- 
separation. The configuration of the coefficient (4.4) is illustrated in Figure 6. 



Figure 6: The rough coefficient a{x) without scale-separation. 


We discretize the spatial domain using a two level mesh as showed in Figure 2, and then we 
solve the optimization problem (3.49) and build the oversampling operator (3.30) on the fine mesh. We 
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Energy Norm 

L 2 Norm 

Numerical solution 

4.44 x W- 2 

1.98 x lO -3 

Corrected numerical solution 

~ 3T7^Tcr^ 

1.18 x llT 3 


Table 2: Numerical errors 


truncate the singular value decomposition of the oversampling operator to e = H. After the offline stage, 
the average number of multi-scale basis functions associated with each edge of the coarse mesh, (4.1), is 
k e ~ 1.00. The smallness of k e reflects the compactness of the local solution space, and the efficiency of 
our method in the offiine stage since the stiffness matrix is small and sparse. 

In the online stage, we choose the forcing function f(x) to be same as (4.2), and measure the error 
of the online numerical solution using the fine mesh solution as reference. The numerical solutions are 
plotted in Figure 7. We can see that the errors in the online numerical solutions are small. We measure 
the error (4.3) in the energy norm and the L 2 norm. The results are summarized in Table 2. Again 
we see that our method achieves very high accuracy in the online stage, which reflects that the good 
performance of our method does not depend on the scale-separation of the coefficient. 



(a) The fine mesh solution. 



(b) Error in the numerical solution, (c) Error corrected using the bubble. 


Figure 7: Online numerical solutions 


4.3 An example with high-contrast channels 

In this subsection, we consider an example with high-contrast channels. The high contrast in the coeffi¬ 
cient violates our uniform ellipticity assumption (1.2), and brings in additional difficulty. The coefficient 
that we consider here is the one with multiple scales (2.14) added with some high conductivity patches 
and channels. log 10 a(x) is plotted in Figure 8, which has very strong heterogeneity. 



Figure 8: log 10 a(ai). The coefficient with high-contrast channels. 

We discretize the problem in the spatial domain as the previous two examples, and build the over- 
sampling operator for each edge e of the coarse mesh. We truncate the singular value decomposition 
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Energy Norm 

L 2 Norm 

Numerical solution 

3.67 x 10“ 2 

1.64 x 10- a 

Corrected numerical solution 

xTTcr^ 

6.13 x 10~ 4 


Table 3: Numerical errors 


of the oversampling operators to e = H. The average number of multi-scale basis functions associated 
with each edge is k e = 0.89 (4.1). Namely, on average, we use less than one multi-scale basis function 
for each edge of the coarse mesh, which implies the compactness of the solution space on local regions 
of the domain for this problem with high contrast coefficient. 

In the online stage, we choose the forcing function to be (4.2), and use fine mesh solution as reference. 
The numerical errors are plotted in Figure 9. We can see that our numerical solutions have high accuracy 
and can capture the large-scale properties of the solution. The numerical errors of the solutions (4.3) 



(a) The fine mesh solution. 


(b) Error in the numerical solution, (c) Error corrected using the bubble. 




Figure 9: Online numerical solutions. 

are listed in Table 3. Again, we see the high accuracy of our online numerical solution. 


5 Concluding Remarks 

In this paper, a novel multi-scale finite-element method is proposed, which is based on the compactness of 
the solution space on local regions of the spatial domain, and does not depend on any scale-separation or 
periodicity assumption of the coefficient. By introducing a set of primary interpolation basis functions, 
we reduce the problem of approximating the local solution space to approximating the trace of the 
solution on each edge of the coarse mesh. We construct basis functions for each edge of the coarse 
mesh separately employing an oversampling operator which is local and compact. The optimal primary 
interpolation basis is also identified by solving a series of under-determined least square problems. 

The resulting method involves two stages: in the offline stage we identify the local compact structure 
of the solution space; in the offline stage we solve the problem efficiently exploiting this compact structure. 
The offline computation is expensive but it only requires solving local problems and is parallel in nature. 
We can rigorously control the error in the online numerical solutions by thresholding in the offline stage. 
The resulting basis functions are called multi-scale basis because they are o(i)-harmonic on each coarse 
element, which have optimal approximation property. Multi-scale basis functions are orthogonal to the 
bubble part of the solution space, and this gives us the flexibility to put back the bubble part of the 
solution in regions where high accuracy is desired. 

Numerical results suggest that our method is very robust and can achieve high accuracy for the 
challenging problems without scale-seperation, or have high-contrast inclusions. 
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Figure 10: The oversampling region. 


APPENDIX 


Proof of Theorem 3.1 


Recall that the restriction operator Pvv->e (3.27), maps u(x)\w to H 1 ' 2 (e) fl C“(e), 

Pw-¥e : Uw(x) ^ u e (x) e H 1/2 (e) rC a (e). (.1) 

And we have introduced the following inner product (3.32) on the domain of Tw~> e , which denote by 
Vw = D(Tw- s-e), 

\\ u w(x)\\v w = a(u] v (x),u] v (x)) + \\uw{x)\\ 2 L 2 {w) + ||div(o(x)Vuw(x))|| z , 2 (M / ) . (.2) 

Note that u e (x) does not necessarily vanish on the two endpoints of e, thus the inner product we defined 
in (3.31) does not apply to v e (x). We first extend u e (x) to the boundary of its two neighborhood coarse 
elements D u , D d , Figure 10. To be specific, we let u e (x) vanish on * i3 , x i4 , x i5 , x i6 , and be piecewise 
linear on d(D u U D d ). Then we solve local boundary value problems on D u and D d separately to extend 
u e (x) to D u U as o(x)-harmonic functions. With this o(a:)-harmonic extension to D u and D d , we 
introduce the following inner product, which agrees with (3.31) for u e {x) that vanishes on xu and Xi 2 . 


||UeO »)" 2 




VM e (i) f a(i:)V« e (s;)di:. 


ID u UD d 


(.3) 


And (.3) can be viewed as an extension of (3.31). 

Next we consider Pi and P2, which are the bounded linear functional that maps uw(x) 6 Vw to its 
values on x 2l and x 22 separately. We can easily see that Pi and P-2 are linearly independent, thus there 
exist <j>i 1 (x ) and cj>i 2 {x), such that 


p Mx ik {x)) = 5j k , j,k= 1,2. (.4) 

We denote the intersection of the kernel of P\(x) and P 2 as Vjy, which is closed subspace of Vw- And 
we denote the projection of ipi^x) and ipi 2 {x) to the orthogonal complement of V^, ( Vw) ± , as ip^ix) 
and ip* (x), namely, 

Ai{x) - ip-^x), i> i2 (x) - ip i2 (x) - ip* 2 (x) _L Vw, (-5) 

where the orthogonality is in the sense of (.2). Then we have 

V’i, (x ih ) = Sj k , \\i’i j {x)\\ Vw >\\il>* j {x)\\v w , j,k = 1,2. (.6) 

And we finish the proof of the existence of unique solutions ^ ( x ) and 1/>* 2 (x) to the optimization problem 
(3.49). Next we proof the property (3.50). 

We choose ’4’* 1 {x)\e and i/>i 2 {x)\ e as the interpolation basis in (3.3), and get the oversampling operator 
Pqs- For any other two interpolation basis functions, i/Uii x ) and ?/>j 2 (a;), we denote the corresponding 
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oversampling operator as Pos, then we compare their singular values, cr)*, and 04 . We use the following 
characterization of singular values, 

CT fe =sup inf ||P OS M(a:)|| H i/ 2 (e) . (.7) 

Vk \\ u ( x )\\v w =l 

where 14 - is a fc-dimensional subspace of Vw- 

For any u(x) £ Vw, we consider its decomposition in V^ and (1 / w) ± , and denote them as Ui(x), 
u 2 (x), 

u(x) = ui(x) + u 2 (x). (.8) 

Then according to our choice the interpolation basis function, we have 

P6s( u l( x )) = P 6sM x )) = Pw->e(u 2 ( x ))- (-9) 

Then according to (.9), we get 

<4 = slip inf ||P 5 s w(z)||tfi/ 2 (e) =sup inf \\P w ^ e u(x)\\ H i/ 2 (e) , (.10) 

14 ll“(*)llv n , = l V k \\H x )\\v w =l 

where 14 - is fc-dimensional subspace of Vfy, since Pq S (( 1 / ^) ± ) = 0. And for at, we have 

a k SU P„ , i? f ||Posw(ar)||Hi/ 2 ( e ) = sup inf \\P w ^eu(x)\\ H i/ 2 {e) . (.11) 

14 \W(x)\\ Vw =i v k ll»(*)llv w -=i 

With (.10) and (.11), we finish the proof of this Theorem. 


Proof of Lemma 3.1 


Denote fi(x) as the piecewise constant functions on the coarse mesh f 1 = L)fL 1 Di, then we consider u,(x ) 
as the corresponding solution to the Poisson equation on fi. Namely, 


-A Ui(x) = fi{x), Ui(x) |an = 0. 


(. 12 ) 


Clear Ui(x) are linearly independent, and we consider the projection of u(x), which is the corresponding 
solution to the Poisson equation with f(x), in the Hg(Q) norm. We denote the projection as Uh(x) = 
CiUi(x), then we have 

(u{x) - u h (x),Ui{x)) H i (fJ) = 0. (.13) 


Using (.13), we have 


IK*) - u h( x ) ||^i (f2 ) = (u{x) - u h (x),u{x) - u h (x)) H i {n] = (u(x),u{x) - u h (x)) H i (n) 

= (u(x) -u h (x),f(x)) L 2 (n) < ||u(®) - u h (x)|| i 2 (f 2 ) ||/(x)|| i 2 (fl) . (.14) 

According to (.13) and (.12), we have 


0=1 (u(x) - Uh(x))fi(x)dx = / (u(x) - Uh{x))dx. 

J O J Di 

Then based on (.15) using the Poincare Friedrichs inequality, and a scaling argument, we have 

[ (u(x) — Uh{x)) 2 dx < C'H 2 f |V(u(®) — Uh(x))| 2 da:.2d:i:. 

J D, J D, 


Summing (.16) over Dj , we have 


||u(x) - u h {x) || L 2 (f2) < CH\\u(x) - u h (x)\\ H i^ n y 
Putting (.17) in (.14), we have 


IK*) - «h(*)llffi(n) ^ CH\\f{x)\\ LH n)\\u(x) - u h {x)\\ H i (n) . 
Then employing ||A(u(x) - «h(a:))||H-i(n) = ||u(*) - u h (x)\\ H i^, we finish the proof 

N 

II fix) -^ c i}i{.x)\\ H -i (n) < CH\\f( x)\\l?(si)- 
2—1 


(.15) 

(.16) 

(.17) 

(.18) 

(.19) 
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